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on the popular Smagorinsky closure. Our results show that, in the context 
of this particular problem, the scale-dependent dynamic modeling approach 
is extremely useful, and reproduces several cstablised results (e.g., the sur- 
face layer similarity theory) with fidelity. Results from both the SGS base 
models are generally in close agreement, although we find a consistent su- 
periority of the Smagorinsky-based SGS model for predicting the inertial 
range scaling of spectra. 

Key words atmospheric boundary layer, large-eddy simulation, neutral, 
subgrid-scale, turbulence. 

Abbreviations ABL - Atmospheric boundary layer; LES - Large-eddy 
simulation; SGS - Subgrid scale; NBL - Neutral boundary layer; TKE - 
Turbulence kinetic energy. 

1 Introduction 

The dynamic subgrid-scale (SGS) modeling approach of Germano et al. 
has been quite successful in large-eddy simulations (LESs) of various engi- 
neering flows 2_ . In this approach, one dynamically computes the values of 
the unknown SGS coefficients at every time and position in the flow. By 
looking at the dynamics of the flow at two different resolved scales, and as- 
suming scale similarity as well as scale invariance of the SGS coefficients, one 
can optimize their values. Thus, the dynamic modeling approach avoids the 
need for a priori specification and tuning of the SGS coefficients. A recent 
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study |31 based on extensive database analysis further suggests that the dy- 
namic modeling approach closely reproduces the minimal simulation error 
strategy (termed as optimal refinement strategy), which is highly desirable 
in turbulence modeling. 

In atmospheric boundary layer (ABL) turbulence, where shear and strat- 
ification and associated flow anisotropics are (almost) ubiquitous, the inher- 
ent scale-invariance assumption of the original dynamic modeling approach 
breaks down. Porte- Agel et al. 0] relaxed this assumption and introduced 
a scale-dependent dynamic modeling approach in which the SGS coeffi- 
cients are assumed to vary as powers of the LES filter width (Af). The 
unknown power-law exponents, and subsequently the SGS coefficients, can 
be determined in a self-consistent manner by filtering at three levels 
In the simulations of neutral boundary layers (NBLs), the scale-dependent 
dynamic SGS model was found to exhibit appropriate dissipation behavior 
and more accurate spectra in comparison to the original (scale-invariant) 
dynamic model Recently the scale-dependent dynamic modeling ap- 

proach was modified and extended by incorporating a localized averaging 
technique in order to simulate intermittent, patchy turbulence in the sta- 
bly stratified flows In parallel, scale-dependent dynamic SGS models 
based on Lagrangian averaging over fluid flow path lines were developed 
by Bou-zeid et al. [8] and Stoll and Porte-Agel 9 to simulate neutrally 
stratified flows over heterogeneous surfaces. 
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The scale-dependent dynamic modeling approach and its variants so far 
always used the popular eddy-viscosity formulation of Smagorinsky |10j as 
the SGS base model. However, this SGS model assumes that the energy dis- 
sipation rate equals the SGS energy production rate. In order to avoid this 
strong assumption, Wong and Lilly 1 1 1 1 proposed a new SGS model based on 
Kolmogorov's scaling hypothesis. A dynamic version of the Wong-Lilly SGS 
model to some extent outperformed the dynamic Smagorinsky model in sim- 
ulations of the buoyancy- driven Rayleigh-Benard convection Further- 
more, the dynamic Wong-Lilly SGS model is computationally inexpensive 
in comparison to the dynamic Smagorinsky SGS model. The combination 
of lesser assumptions and cheaper computational cost certainly make the 
Wong-Lilly model an attractive SGS base model for LES. Therefore it is of 
interest to explore if the Wong-Lilly SGS model or its variants are capa- 
ble of simulating different flow regimes of the ABL. It is generally agreed 
upon that in comparison to buoyancy-driven flows, large-eddy simulations 
of shear-driven boundary layer flows are far more challenging. Thus, in the 
present study, we focus on neutrally buoyant shear-driven ABL flow. In 
order to realistically account for the near-wall shear effects, we first formu- 
late a locally averaged scale-dependent dynamic version of the Wong-Lilly 
SGS model (henceforth LASDD-WL, see Appendix for details). Then, we 
comprehensively compare its performance with the locally averaged scale- 
dependent dynamic Smagorinsky (hereafter LASDD-SM) SGS model earlier 
developed by Basu and Porte- Agel UJ. 
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The structure of this paper is as follows. In Section 2, we briefly provide 
the technical details of a case study. Extensive comparisons (in terms of the 
similarity theory, spectra, and flow visualizations) between the LASDD-WL 
and LASDD-SM SGS models arc performed in Section 3. Finally, concluding 
remarks are made in Section 4. 

2 Description of Simulations 

In this work, we perform large-eddy simulations of a turbulent Ekman layer 
(i.e., pure shear flow with a neutrally stratified environment in a rotating 
system) utilizing the LASDD-SM [H10 and LASDD-WL (see Appendix) 
SGS models. Both these simulations are identical in terms of initial con- 
ditions, forcings, and numerical specifications (e.g., time integration, grid 
spacing). Technical details of our LES code and the LASDD-SM SGS mod- 
eling approach have been described in detail in jo] and will not be repeated 
here for brevity. 

The selected case study is similar to that of the LES intercomparison 
study by Andren et al. |12| . The simulated boundary layer is driven by an 
imposed geostrophic wind of (U g , V g ) = (10, 0) ms _1 . The Coriolis parame- 
ter is equal to f c — 10 -4 s _1 , corresponding to latitude 45° N. The compu- 
tational domain size is: L x = L y = 4000 m and L z — 1500 m. This domain 
is divided into N x x N v x N z = 40 x 40 x 40 nodes (i.e., A x = A y = 100 
m, and A z = 38.5 m). The motivation behind the selection of this coarse 
grid-resolution is two-fold. Primarily it allows us to perform a direct compar- 
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ison with the results from 12 , which used almost the same grid-resolutions. 
More importantly, coarse grid-resolution enables us to identify the strengths 
and/or weaknesses of different SGS models, as well as, to underscore their 
impacts on large-eddy simulations. The simulations are run for a period of 
10 x f- 1 (i.e., 100,000 s), with time steps of 2 s. The last 3 x /- 1 interval 
is used to compute statistics. A passive scalar is introduced in the flows by 
imposing a constant flux (wcq) of 10 -3 kg m~ 2 s -1 at the surface. The lower 
boundary condition is based on the Monin-Obukhov similarity theory with 
a surface roughness length of z =0.1 m. 



3 Results and Discussions 

In this section, we report the results of the LASDD-SM and LASDD-WL 
SGS models-based simulations and compare them with results from the in- 
tercomparison study |12j . wherever possible. This particular case (without 
the inclusion of passive scalars) was also simulated by Kosovic |13| using a 
nonlinear SGS model, and recently by Chow et al. (T3|, who utilized a so- 
phisticated hybrid SGS model. Our simulations show that both the LASDD 
SGS models perform very well, and the results are comparable to the past 
studies. 

Temporal evolution of the surface friction velocity (it*) is very similar in 
both the simulations (not shown). The average value of u* during the last 
3 x interval is approximately 0.44 ms" 1 in the case of the LASDD-SM 
model. The LASDD-WL model produces a marginally higher value (0.454 
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ms -1 ). The corresponding values found in ^2] are: 0.425 ms _1 (Moeng), 
0.448 ms -1 (Mason - backscatter), 0.402 ms _1 (Mason - non-backscatter), 
0.402 ms" 1 (Nieuwstadt), and 0.425 ms _1 (Schumann). 
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Fig. 1 Temporal evolution of the nonstationarity parameters C u (left) and C v 
(right). 



In Figure ^ we present the nonstationary parameters C u and C v (see 
[T2^ for definitions). Under steady state conditions, these parameters should 
approach unity. Although none of the past |12II14| and present simulations 
are quite close to steady state conditions, they are more or less in phase 
with each other. All these simulations clearly portray the inertial oscillation 
of period 2tt// c , as anticipated. 

Accurately simulating the non-dimensional velocity gradient (4>m), and 
the scalar gradient (4>c) rn the neutrally stratified surface layer has proven 
to be a very challenging task for many atmospheric LES models. It is well 
known that the traditional Smagorinsky model is over-dissipative in the 
near-surface region and gives rise to excessive mean gradients in velocity 
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and scalar fields (cf. 52)- Fortunately, state-of-the-art LES-SGS modeling 
approaches of Mason and Thomson Sullivan et al. ^Sjj Kosovic |LS| . 
Porte- Agel et al. |U, Porte- Agel 0, Esau [T7], Chow et al. [Hj, Bou-zeid 
et al. [Hj, and Stoll and Porte- Agel 9 offer major improvements over tradi- 
tional Smagorinsky-type SGS models, and reproduce the non-dimensional 
gradients reasonably well. From Figure it is clear that both the LASDD- 
SM and LASDD-WL SGS models behave satisfactorily, albeit, the perfor- 
mance of the LASDD-WL SGS model is superior. We would like to stress 
that both the LASDD SGS modeling approaches do not require any addi- 
tional stochastic term or supplementary near-wall stress models for reliable 
performance in an LES. In the framework of Monin-Obukhov similarity the- 
ory, the non-dimensional velocity gradient (</>m) is indisputably equal to one 
(the dotted line in Figure 0- left). However, in the literature there is no 
consensus on the 'true' magnitude of the non-dimensional scalar gradient 
(<fic)- Businger et al. ^H|> based on the Kansas field experiment, proposed 
a value of 0.74. Recent field observations, however, suggest values close to 
0.9 (for a review, see |19|L From the present coarse-resolution simulations, 
it is difficult to favor either of these values. However, qualitatively, both the 
LASDD SGS models portray very similar non-dimensional scalar gradient 
profiles (Figure 0- right). 

In neutrally stratified ABL flows, the observed peak normalized velocity 
variances occur near the surface and are of the magnitude: ct^/m^ ~ 5 — 7, 
a vl u 1 ~ 3 — 4, and a^/ul ~ 1 — 2 [201 - The corresponding values found 
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Fig. 2 Simulated non-dimensional velocity (left) and scalar (right) gradients. The 
dashed lines correspond to the values of 1 (left) and 0.9 (right). These values are 
expected to hold in the surface layer under neutral conditions according to the 
similarity theory. 

in our simulations (Figure approximately fall in these ranges. The sim- 
ulated results also concur with the outer boundary layer observations. For 
example, the KONTUR data |2U give al/ul ~ vl/ut ~ 1 and a* /v* ~ 0.5 
at z = 0.75^ (where Zi denotes the inversion height). The normalized scalar 
variances (c;?/c*) are also shown in Figure OH Here, c* is the surface scalar 
scale (— — wcq/u*). In |12j . it was found that the consensus among dif- 
ferent SGS models is poorer in the case of passive scalar in comparison 
to the momentum case. The disagreements between different SGS models 
could be partially attributed to different a priori prescriptions for the SGS 
Prandtl (Prscs) number, and underscore the need for the determination 
of Prscs in a self-consistent manner, as is done in the present study. One 
must also acknowledge the facts that the passive scalars exhibit complex 
spatio-temporal structure, and the statistical and dynamical properties of 



10 



Anderson et al. 



passive scalars are remarkably different from the underlying velocity fields 
[52123] . 

We point out that the individual plots in Figure |3J represent both the 
normalized resolved and total (resolved + SGS) variances. In the LASDD 
modeling approach, one does not solve additional prognostic equations for 
the SGS turbulence kinetic energy (TKE) and the SGS scalar variances. 
However, the SGS variances can be roughly diagnosed using the approach 
of Mason [23] . 
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Fig. 3 Simulated normalized longitudinal (top-left), transverse (top-right), ver- 
tical (bottom-left) velocity variances. Simulated normalized scalar variances are 
shown in the bottom-right plot. 
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The one-dimensional longitudinal velocity and passive scalar spectra 
are computed at heights z = O.lzi, and z = 0.5zi, and presented in Fig- 
ure 21 The spectra highlight the most important difference between the 
LASDD-SM and LASDD-WL SGS models: the LASDD-WL SGS model 
seems to be over-dissipative (indicated by steeper spectral slopes at higher 
wavenumbers). In the case of the LASDD-SM model, the longitudinal veloc- 
ity and scalar spectra clearly show extended inertial range (k x 5 ^ 3 scaling) 
at z = 0.5zi . Near the surface {z = O.lzj), the longitudinal velocity spec- 
tra show the anticipated production range (fcj -1 ), as well as a short inertial 
range. Recent research suggests that the production range is (likely) related 
to elongated streaky velocity structures (see below). Traditional SGS models 
typically do not reproduce well defined inertial ranges in coarse-resolution 
simulations (cf. 12 ). From that perspective, the performance of the LASDD 
models could be considered a near success. We note that the original plane- 
averaged 0] and the Lagrangian-averaged 8,9^ scale-dependent dynamic 
SGS models also reproduced the characteristics of the one-dimensional lon- 
gitudinal velocity spectra remarkably well. However, near the surface, the 
passive scalar spectra predicted by these SGS models showed unphysical 
pile up of scalar variances (50 . This was possibly due to small dynamically 
determined eddy-diffusion coefficients near the surface 5 9. In the present 
study we did not encounter this issue. 

A few previous LES studies have reported the existence of elongated 
streaky structures in the neutral surface layers ^3,16,25,26 27,28 . The 
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Fig. 4 Spectra of longitudinal velocity (top), and passive scalar (bottom). The 
spectra are plotted for z = O.lzi (left), and z = 0.5zi (right) levels. The dashed 
and dotted lines depict the inertial range (fc^~ 5 ^ 3 ) and production range (fc^ 1 ) 
scalings, respectively. 



link between experimentally observed long production range (fc _1 scaling) 
in the streamwise spectra of the longitudinal velocity and the elongated 
streaky structures has recently been discussed in depth by Carlotti |28|. 
Moreover, strong correlations between these streaky structures and large 
negative momentum flux were earlier reported by |26j . From Figure [3] (top), 
it is clear that both the LASDD models show streaky structures, roughly 
parallel to the mean wind direction, in the surface layer (at z = O.lZj). 
However, significant morphologic differences are noticeable in the mid-ABL 
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flow structures. In accordance with past studies (cf. j2E])i the LASDD- 
SM SGS model predicts non-coherent structures at z = 0.5z,i. In contrast, 
large coherent structures persist in the LASDD-WL model results (Figure|5l 
bottom-right). Another interesting feature of this plot is the (virtual) non- 
existence of fine-scale flow structures. This can be directly associated with 
the over-dissipative nature of the LASDD-WL SGS model, as discussed be- 
fore. In essence, we can infer that the (non-)existence of coherent structures 
in NBL flows are strongly dependent on SGS parameterizations, especially 
for coarse-resolution simulations. A few previous studies somewhat support 
this inference. For instance, the nonlinear SGS model |Li| . and the mod- 
ified Smagorinsky SGS model [27J barely produced any elongated streaky 
structures. 

4 Concluding Remarks 

Two locally averaged scale-dependent dynamic SGS closures - the LASDD- 
SM jHHZj and the LASDD-WL (this study) - have been used to simulate a 
neutral ABL case. Although the theoretical foundations of these SGS mod- 
els are fundamentally different, results presented in Figures ^ through El 
illustrate strong congruence between their results, and with firmly estab- 
lished results (i.e. the Monin-Obukhov similarity theory and the inertial 
range scaling of spectra). The normalized variances computed in our simu- 
lations also closely follow the ones calculated from field measurements. The 
major noticeable and consistent difference between the results is shown in 
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Fig. 5 Visualization of longitudinal velocity fields simulated by LASDD-SM 
(left), and LASDD-WL (right) SGS models. The horizontal cross-sections are 
taken at z = O.lzi (top), and z — 0.5zi (bottom). 

Figure the LASDD-WL SGS model appears to be over-dissipative at the 
higher wavenumbers, in comparison to the LASDD-SM SGS model. In Fig- 
ure we see that both SGS models predict elongated streaky structures 
in the near-wall region (z = O.lzi). These coherent structures are no longer 
evident at higher locations in the domain, in the case of the LASDD-SM 
SGS model-based simulation. Due to undue dissipation, the LASDD-WL 
SGS model-based simulation results in quite different flow structures at 
this level. 
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The Wong-Lilly SGS base model requires fewer assumptions and comes 
at slightly less computational cost in comparison to the commonly used 
Smagorinsky SGS base model. Unfortunately, these advantages seem to be 
offset by its over-dissipativc tendency at higher wavenumbers. Some inher- 
ent assumptions of the Smagorinsky base model can also be eliminated by 
solving a prognostic equation for the TKE. However, when using this TKE 
SGS approach, the SGS model coefficients are often tuned for different ABL 
flow conditions |16U29| . An alternative approach would be to formulate a 
dynamic version of the TKE SGS model, which will also account for en- 
ergy backscatter. We are currently working on this SGS approach to better 
represent the physics of atmospheric boundary layer flows. 

Acknowledgements This work was partially funded by the National Institute 
of Standards and Technology, the National Science Foundation and the Texas 
Advanced Research Program grants. All the computational resources were kindly 
provided by the High Performance Computing Center at Texas Tech University. 

Appendix 

The SGS model proposed by Wong and Lilly |llj can be written as: 

Tij — -TkkSij = — 2CwL^J i Sij (Al) 
6 ' 

where Tij and Sij denote the SGS stress tensor and the resolved strain rate tensor, 
respectively. Cwl is a model coefficient to be specified or determined dynamically. 
In a recent LES study of neutral boundary layer flows, Chow et al. |14| utilized 
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a dynamic version of this SGS model in conjunction with the approximate de- 
convolution model (ADM) for resolvable subfilter-scale (RSFS) components. To 
account for the smaller underresolved eddies in the surface layer, they used a near- 
wall stress model in addition to the dynamic Wong-Lilly SGS and ADM-RSFS 
models. As an alternative approach, in this work, we formulate a locally averaged 
scale-dependent dynamic version of Equation (named LASDD-WL). 

The SGS stress tensor (Tij) at the filter scale (Af) is defined as: Tij = 
UiUj — UiUj. In a seminal work, Germano et al. pQ proposed to invoke an ad- 
ditional explicit test filter of width aAf in order to dynamically compute the 
SGS coefficients. Consecutive filtering at scales Af and at aAf leads to a SGS 
turbulent stress tensor (Tij) at the test filter scale aAf. 



Tij = m uj — Ui Uj, (A2) 



where an overline (■ ■ ■) denotes filtering at a scale of aAf. From the definitions 
of nj and Tij an algebraic relation can be formed, known in the literature as the 
Germano identity: 



Lij — UiUj Ui Uj — Tij Tij. (A3) 

This identity is then effectively used to dynamically obtain unknown SGS model 
coefficients. In the case of the Wong-Lilly model (Equation HAlll 'l. this identity 
yields: 

Lij - ^L kk Sij = {Cwh) Af Mij, (A4) 
4 / 3 -.4/3 < - CwL) °'*f 



2Af ^ 1 — a (C W l) — J ^ one assumes sca l e invariance, 
i.e., {CwL) aAf = (CwL) Af , then the unknown coefficient [CwL) Aj can be easily 
determined following the error minimization approach of Lilly |3U| : 
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In the context of the present study, the angular brackets {• ■ ■) denote localized 
spatial averaging on horizontal planes with a stencil of three by three grid points 

m 

Recent studies have shown that the assumption of scale invariance is seriously 
flawed for sheared and stratified boundary layer flows | -115161 7 8 9 . In other words, 
the ratio of (Cwh) aAf to (Cwl) A{ should not be assumed equal to one for 
most of these ABL flow scenarios. Rather, this scale-dependence ratio should 
be determined dynamically. In order to implement the scale-dependent dynamic 
procedure, one needs to employ a second test filtering operation at a scale of a 2 Af 
[denoted by (■■■)]. Invoking the Germano identity for the second time leads to: 

Qij - -QkkSij = (C W L) Af N ij , (A6) 



where 



Qij — UiUj UiUj 



and 



This results in: 



Ni > =2A * y- a (c WL ) Af ) **• 

^=§Sf (A7) 

Following |U, the following scale-dependence assumption can be made: 

{CwL) aA (CwL) a 2 A 

P = (r ^ = (r \ L > ( A8 ) 

This is a much weaker assumption than the scale-invariance modeling assumption 
of (3 — 1. Now, from Equations l|A50 and l|A7^ . using Equation l|A8fl . one solves 
for the unknown parameter j3, which in turn is used to compute the Wong-Lilly 
SGS model coefficient, (CwL) Af , utilizing Equation 1A5I . 
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Solving for ft essentially involves finding the roots of a fifth-order polynomial 

4: 

A + Ai/3 + A 2 (3 2 + A z fi A + A 4 (3 4 + A 5 (3 5 = (A9) 

where Ay = 0403 — a^og, Ai = a\04 — a-jag, A2 = 0203 + aias — a^ag, A3 = 
0,20,4, — 0709, A4 = 0205 — agaio, and As = —07010. In the case of Wong-Lilly 
SGS base model, we derive: a\ = {QijSij), 02 = (— a 8 ^ 3 QijSij), 03 — (Sij ), 04 = 
(-2a 4 / 3 ^- 2 ), a 5 = (a 8 / 3 ^ 2 ), a 6 = (L«S£), a 7 = (-a 4 / 3 !^), a 8 = }, 
ag = (~2a 8 ^ 3 Sij }, and aio = {a 16//3 5'i :) }. Please note that the coefficients (ai to 
aio) involve significantly lesser number of tensor terms in comparison to the ones 
derived by Porte-Agel et al. [I] using the Smagorinsky SGS base model. Lesser 
number of calculations (specifically the tensor multiplications) undoubtedly lead 
to cheaper computational costs. 

Scale-dependent formulation for scalars can be derived in a similar manner 
5. . The Wong-Lilly model for a generic scalar (c) could be written as: 

Cwl .4/3 dc . . , 

* = -P^ A * fl* (A10) 

where Ptsgs is the so-called SGS Prandtl number. In the dynamic or scale- 
dependent dynamic modeling approaches, typically the lumped SGS coefficient 
(CwLPr],Q S ) is determined in a self-consistent manner. This procedure not only 
eliminates the need for any ad hoc assumption about the SGS Prandtl number 
{PrsGs), h also completely decouples the SGS scalar flux estimation from SGS 
stress computation. In the scale-dependent approach one further defines a 
scale-dependent parameter for scalars analogous to Equation l|A8fl . For the 

Wong-Lilly SGS base model, it could be written as: 

(C WL Pr-^ s ) aA (CwLPr- 1 GS ) a2A 

ft = -7 =1 r = -7 rr^ — -i ( An ) 

{Cw L Pr SGS ) {CwLPr SGS ) aA 
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As before, f3 c could be determined by solving the fifth-order polynomial: 

A Q + A^c + A 2 (f c + Ai0i + AiPi + A 5 f3l = (A12) 

where Aq = 0103 — aeag, A\ = ai<i4 — 070%, A2 = 0203 + oias — a^ag, A3 = 
12^4 — 0,70,9, A4 — 020s — agOio, and A$ = — 0,70*10. For the Wong-Lilly SGS 
base model for scalars, we get: 01 = {K^-g^-), 02 = (—a^^K'i^-), 03 = }, 

<* = (~^ 4/3 S- 2 )' «s = <« 8/3 £ 2 ), o 6 = (K t §-), a 7 = (-a^Kify, o 8 = 
2 — ^ 2 2 

(^7 >. fl 9 = <- 2a8/3 ^7 >> and 010 = (a 16/3 ^7 >• Here, = (u^-Z- fj, and 
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